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ABSTRACT 

Convenient and efficient function generation subprograms have 
been developed for handling functions of one variable and two types 
of functions of two variables. These subprograms can easily be used 
In any digital or hybrid simulation requiring function generation. 
Use of these programs can often lower overall program execution 
time. 
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FUNCTION GENERATION SUBPROGRAMS FOR USE IN DIGITAL SIMULATIONS 
By Cl int E. Hart 
Lewis Research Center 

SUMMARY 

Dynamic models of propulsion systems often contain component 
performance data in the form of functions of one and two variables. 
Digital and hybrid simulation of these dynamic models require sub- 
programs to handle function generation. Convenient and efficient 
function generation subprograms have been developed for handling 
functions of one variable and two types of functions of two 
variables. Use of these subprograms in digital or hybrid simula- 
tions can often lower the overall program execution time. 
INTRODUCTION 

Simulation of dynamic models of turbojet or tubofan engines 
has proved valuable in analyzing their dynamic behavior and in 
designing engine control systems. Often the model equations con- 
tain Several functions of one and two variables* The two variable 
functions are difficult to program accurately for analog computer 
simulation and require large amounts of analog ccxnputer equipment. 
For these reasons most engine simulations have been programmed for 
digital computers. 

The most commonly used programming language for these engine 
simulations Is FORTRAN. However, some simulations have been 
programmed in either CSMP or one of the special languages developed 



for hybrid computers. When CSMP first became available It had no 
provision for handling a function of two variables. Also, at that 
time no suitable function generation subprograms were being retained 
in our computer system library. Thus, to use CSMP for engine simula- 
tions, two variable function generation subprograms had to be 
developed. 

For FORTRAN engine simulations the programmers usually wrote 
their own subprograms to handle functions of one and two variables. 
Examination of many of these subprograms showed 

that they were much slower in execution time than desirable. Fast 
execution times for these subprograms is essential because in some 
dynamic engine simulations the function generation subprograms may 
be called more than 10^ times per second of problem time. 

Therefore, a progranwing effort was made to develop convenient 
and efficient one and two variable function generation subprograms. 

The following sections will contain discussion of the types of 
functions usually encountered and descriptions of the function 
generation subprograms. The application of these function genera- 
tion subprograms with CSMP or other digital simulation programs will 
be explained. 

PROGRAM DEVELOPMENT 

A typical function frequently encountered in engine simulations 
is a function of one variable Z = f(X), as shown In figure I. For 
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the simulation, a set of discrete points Is selected to represent this 
function. Usually linear interpolation is used between points, but 
higher order Interpolation can be used if desired. The X-input values 
are often called breakpoints when using linear interpolation. 

There are two types of functions of two variables that may occur 
in engine simulations, A common type of function is shown in figure 
2. For this type the range of the family of curves Is over the same 
set of values of the X-input variable. Usually, two dimensional linear 
interpolation is made in directions parallel to the X-input and Z-outnut 
axes, 

A second type of function of two variables is shown in figure 3. 
This type, often called a map, differs from the first type in that 
the range of each curve is over a different set of values of the 
X-inout variable. The two-dimensional linear I nterpola t Ion ment ioned 
above cannot be used for this map-type function of two variables. 
Significant errors can occur with this interpolation method, A 
special method should be used with interpolation in directions not 
parallel to either the X-input or Z-output axes. An interpolation 
method of this type will be discussed in a later section. 

Description of Subprograms 

Consider first a subprogram to handle a function of two variables 
as shown in figure 2. The data which describes the function consists 
of lists or tables of X breakpoints and Y curve values and the 
correspond i ng Z output values. For this type of function the X 
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breakpoint values are common for all curves. The X and Y values are 
listed in order of monoton i cal ly increasing value. The Z values are 
listed in order corresponding to the X values for the first curve 
(lowest Y value), second eurve, etCt 

Typically, a comparison table search method is used in function 
generation subprograms. The X and Y inputs are compared with table 
values until the proper pairs of table values are found which bracket 
the inputs. Usually the table search starts at either the low or 
high end of the table and moves either up or down. 

The table search method can be made more efficient, i.e., faster, 
by using table-search indices. These search indices are used as start 
ing points in the table searches of the X and Y data tables. First, 
the X and Y inputs are tested to see if they are In the same intervals 
between pairs of X or Y values that they were in for the preceding 
entry. If not, the search indices are incremented either up or do’wn 
until the proper intervals are found. The final values of the search 
indices are stored for use at the next entry. After the proper table 
values are determined a two-dimensional linear Interpolation Is made 
to obtain the output. 

The subprogram for a function of one variable is similar to the 
one described above. Since there is only one curve, no Y values are 
required, and only one set of Z values. A table search index and 
one dimensional linear interpolation are also used. 



The subprogram for a map* type function of two variables requires 
a set of X values for each curve. The special interpolation method, 
which will be described later, requires that the same number of break- 
points must be used to define each curve. The X values are listed 
in order of monoton lea 1 ly Increasing value for the first curve (lowest 
Y value), second curve, etc. The Z values are listed in order 
corresponding to the X values* Additional data required for all of 
the function generation programs are the number of points per curve 
and in the case of the two-variable subprograms, the number of 
cu rves. 

Part of the map-type function of two variables of figure 3 is 
also shown tn figure 4. Extra lines and points have been added to 
explain the Interpolation method which is used in the subprogram. 

Consider the point A whose location is determined by inputs 
and Y^. First, the input Y/\ is compared with successive 
values in the Y table to determine which pair of curves it is 
between. In this example, it is between curves having values of 
Y(2) and Y(3). Adjacent pairs of similarly located breakpoints on 
these curves define quadrilaterals. Through special logic in the sub- 
program the quadrilateral containing point A can be determined. In 
this example it is defined by points B, C, D, and E. Knowing the Z 
and X coordinates of these points, the Z and X coordinates of 
points F and G can be found by linear interpolation in Y along 
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lines BC and ED, Then the Z output 2JS can be found by linear 
interpolation in X along line FG. 

Another problem which must be considered in developing these sub- 
programs is; what to do when one of the inputs Is beyond the range 
of the table data. If an input Is out of range of the tabulated 
values for the function, the output is calculated by extrapolating 
beyond the end values of the tables. An error message is printed, 
but only the first time an input goes out of range. The decision 
whether or not to extrapolate or print an error message is arbitrary 
and the subprograms can easily be modified to fit needs of a partic- 
ular s imulat ion. 

FORTRAN listings for two subprograms, FUN2 and MAPFUN, to 
handle both types of functions of two variables are presented in 
Appendix A. Also in Appendix A is a listing of a subprogram FUNi to 
handle a function of one variable. 

Instructions for Using Subprograms 

The function data should be stored in one-dimensional arrays. 
Separate arrays should be used for each set of X Input, Y input, 
and Z output data. For FORTRAN simulations the BLOCK DATA subprogram 
is convenient to use for entering the function data. Duplicate 
labelled COMMON statements must be used In the main simulation program 
and the BLOCK DATA subprogram. For CSMP simulations the TABLF data 
Input option is used to enter the function data. 



In addition to arrays for function data, two arrays are required 
for the search indices. One of these arrays is used to store the 
X-input search indices and one to store the Y-input search indices 
■for all the functions. Also, an array is needed to store an error 
indicator which controls printing an error message when an input 
Is out of range of the function data. For FORTRAN simulations these 
arrays can be initialized by DATA statements. For CSMP simulations 
a subprogram FUNSET can be called in the INITIAL section to initialize 
these arrays. A listing of this program is presented In Appendix A. 

The calling statement for a function of one variable is of the 

form: 

ZOUT = FUNl (N,NXP,X1 ,Z1 ,XIN) 

N is an Integer used to select the table search index for each 
call. NXP is the number of points per curve, XI and Zl are 
arrays containing the function data, XIN is the input variable and 
ZOUT is the output variable. 

For the regular type function of two variables the calling state- 
ment Is of the form; 

ZOUT = FUN2(N,NXP,NYC,XX1,ZZ1 ,X1N,YIN). 

ZOUT, N, NXP, and XIN are the same as described in the preceding 
paragraph. NYC is the number of curves. XX1,YY1, and ZZ 1 are arrays 
containing the function data and YIN is the second input variable. 

The calling statement for the map- type function of two variables 
is of the form: 


ZOUT 


MAPFUN(N,NXP,NYC,XM1,YM1,ZM1,XIN,YIM). 
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XMl, YMl, ZMl are arrays containJng the function data. All the 
other arguments are the same as described for the regular function of 
two variables. 

CONCLUDING REMARKS 

For digital simulations requiring function generation, considerable 
savings In execution time can be made by using the subprograms discussed 
In this memorandum. Timing tests Indicate these subprograms run at 
least 80% faster than comparable subprograms used In current digital 
engine simulations. These subprograms can be easily Incorporated in 
digital simulation programs. 



APPENDIX - FORTRAN SUBPROGRAM LISTINGS 


FUNCTION Ftmi(N,NXP,XX,ZZ,XIN) 

COMMON /FHFMR/ I X ( 4 0 ) ^ UY( 4 0 ) , 1 PRR C 40 ) 

DIMENSION XX(NXP),ZZ(MXP) 

I = IX(N) 

C*****TEST for X IN PREVIOl'S I NTPR^'^L***** 
IF(XIM-XX(D) 120,200,110 
110 IF(XIN-XX(I+D) 200, 140,140 
C*****COUNT DOV/N***** 

120 IF(XIN-XXd)) 160,160, 130 
130 I = (-1 

IF(XIM-XX(D) 150,200,200 
C** * * *C01J NT UP*i>r*** 

140 IF(XIN-XX(NXP)) 150,170,170 
150 I = 1+1 

IF(XfN-XX(l + D) 200, 200, 150 
160 I « 1 

no TO 1«0 
170 I = NXP-1 

180 IF(IFRR(N)) 200,100,100 
100 WRITE(6,400) M,XIN 
IFRR(N) = -1 

C*++**INTERPOLATE FOR ANSWER***** 

200 XFRAC = (XIN-XXCI ))/(XX( !+1)-XX( m 
FUNl = ZZ(I)+XFRAC*(Z7(I+1)-Z7CI)) 

IX(N) = I 
RETURN 

400 FORMATCIHO, 12MF!iNPTinN NO.,I^/20M ITPUT HUT or 
12X,6HXIN = ,ni2.4) 

END 


FUNCTION FUN2(N,NXP,NYC,XX, YY,ZZ,X IN, YIN) 
COMMON /FMEMR/ I X ( 4 0 ) , J Y ( 40 ) , I FRR( 40 ) 
DIMENSION XX(NXP), YY(NYC),ZZ(NXP,NYC) 

I * IX(N) 

J « JY(N) 

C*****TEST FOR X IN PREVIOUS I NTERVAI ***** 
IF(XIN-XX(D) 120,200, 110 
110 IF(XIN-XX(I + D) 200, 140, 140 
C*****COUNT DOWN***** 

120 IF(XIN-XXd)) 160, 160, 150 
130 I » 1-1 

IF(XIN-XX(I)) 130,200,200 
C*****COUNT UP***** 

140 IF(XIN-XX(NXP)) 150,170,170 
150 I * 1+1 

IF(X|M-XX(I + D) 200,200, 150 
160 I » 1 

no TO 180 
170 I = NXP-1 
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180 IFdRRR(M)) 200,190, 190 
190 V/RITE(6,400) N,XIN,YIM 
lERR(N) = -1 

C****«TEST FOR Y IM PREVIOUS INTERVAL***** 

200 IF(YIN-YY(J)> 220,300,210 
210 IF(YIN-YY(J+1) ) 300,240,240 
C*****COUNT DOWN***** 

220 IF(YIN-YY(D) 260, 260, 230 
230 J = J-1 

IF(YIN-YY(J)) 230,300,300 
C*****OOUNT IIP***** 

240 IF(YIN-YY(NYO)) 250,270,270 
250 J = J+1 

IF(YIN-YY(J + D) 300,300, 250 
260 .1 = 1 

GO TO 280 
270 .) = NYC-1 

280 IFdERR(N)) 300, 290, 290 
290 WRITE(6,400) N,XIN,YIN 
lERR(N) = -1 

c ***** interpolate for answer***** 

300 XFRAC = (XIM-XXd ))/(XXd+l)-XX(n) 

PIZZ = ZZd,J)+XFRAC*(ZZd + l,J)-ZZd,J)) 

P2ZZ - ZZd , J + l)+XFRAC*(ZZd + l, J*l)-ZZd ,.l + l) ) 

YFRAC = (YIM-YY(J) )/(YY(J+l)-YY(d) ) 

FUN2 • P1ZZ+YFRAC*(R2ZZ-P1ZZ) 

IX(N) = I 
JY(N) = J 
RETURN 

400 FORMAT(1HO,12HFUNCTION NO.,I3,20H INPUTS OUT OF RUiNCF, 
12X,6HXIN = ,012.4, 2X,6HYIN = ,012.4) 

END 


REAL FUNCTION MAPFUN CM, NXP, MYC, XVALS, YVALS, ZVALS, V | >1, Y I M ) 
COMMON /FMEMR/ I X (40 ) , J Y( 4 0 ) , I pRRC 40 ) 

DIMENSION XVALS(NXP,MYC), YVALS(NYC),ZVALrCNX'’,NYC) 

I = IXCN) 

J = JY(M) 

C*****TEST FOR Y IN PREVIOUS I mtER'/AI.***** 

IFCYIN-YVALSCJ)) 120,200,110 
110 IFCYIN-YVALSCJ+D) 200, 140,140 

d*****count down***** 

120 IFCYIN-YVALSCD) 160,160, 130 
130 J = J-1 

IFCYIN-YVALSCJ)) 130,200,200 
C*****COUNT IIP***** 

140 I F( YIN-YVALSCNYC) ) 150,170,170 
150 J = J+1 

IFCYIN-YVALSCJ + D) 200,150,150 
160 J » 1 

GO TO 180 
170 J = NYC-1 
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180 IF(IERR(N))2D0,1'!0,190 
190 WRITE(B,400)N,XIN,YIN 
IF.RROJ) = -1 

C*»***nAI.CULATE XLO AND XHI***»* 

200 YFRAC => (YIN-YVALR(J))/(YVAt.S(J + l)-YVALS(>l)) 

Xl.O = XVALS(l,J)+YFRAf;*(XVAI.S(l,J + l)-XVALR(l,J)) 

XHI = XVALS(I+1,J)+VFRAC*(XVALS(I+1,J + 1)-XVA1.R( l + l,.D) 

c*****tf.rt for x between XI.0 ano xhi****» 

IF(XIM-XLO) 220,300,210 
210 IF(XIN-XHI) 300,240,240 
c»»**»r:oiiNT nowN***** 

220 XMIN = XVALSd, J)+YFRA0*(XVAI,R(1, J + 1)-XVAI.R(1,J) ) 
IF(XIN-XMIN) 2B0,260,230 
230 I = 1-1 
XHI = XLO 

XLO = XVALS(l,il)+YFRAr'»(XVAI.S(l,,l + l)-XVAI.R(l,.l)) 
IF(XIM-XLO) 230,300,300 
C*****COUNT UP***** 

240 XHAX = XVALS(NXP,J)+YFRAC*(XVAI.?^(NXP,vl + l)-XVALr>(NXP,J)) 
IF(XIN-XMAX) 250^270,270 
250 I - I+l 
XLO ~ XH I 

XHI = XVALSC 1+1, J)+YFRAC*(XVAI.S( I+1,0 + 1)-XVAL':( l+l,vl) ) 
IF(XIM-XHI) 300,250,250 
260 I = 1 

on TO 280 
270 I = NXP-1 

280 IF(IERR(N))300,290,290 
290 WRITE(6,400)N,XIN,YIN 
lERR(N) « -1 

INTERPOLATE FOR ANSWER***** 

300 XFRAO = (XIM-XLO)/(XHI-XLO) 

7.1. * 7VALS( I ,il)+YFRAO*(2VAI.S(l ,>! + l)-ZVALS( I ,.0 ) 

7R = ZVAI.S( 1+1, J)+YFRAO*(ZVAI.S{ 1 + 1, J + 1)-7.VAI.SC l+l,.l) ) 

MAPFUN = ZL+XFRAC*(ZR-ZI.) 

IX(N) = I 
vlY(N) = J 
RETURN 

400 FORMATCIHO, 8H HAP NO.,I3,20H INPUTS OUT OF RA'IOE, 
12X,6HXIM = ,0,12 .4,2X,6HV|N = ,012.4) 

ENO 


FUNOTIOM FIINSET(N) 

COMMON /FMEMR/ I X (40) , ,!Y( 4 0 ) , I PRR( 40 ) 
00 10 K=1,N 
IX(K) = 1 
JY(K) = 1 
lERR(K) = 1 
10 CONTINUE 

FlINSET = 1.0 

RETURN 

END 
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